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Abstract. We report Molecular Dynamics simulations for a new model of tetrahedral 
network glass-former, based on short-range, spherical potentials. Despite the simplicity 
of the forccfield employed, our model reproduces some essential physical properties of 
silica, an archetypal network-forming material. Structural and dynamical properties, 
including dynamic heterogeneities and the nature of local rearrangements, arc 
investigated in detail and a direct comparison with models of close-packed, fragile 
glass-formers is performed. The outcome of this comparison is rationalized in terms 
of the properties of the Potential Energy Surface, focusing on the unstable modes 
of the stationary points. Our results indicate that the weak degree of dynamic 
heterogeneity observed in network glass-formers may be attributed to an excess of 
localized unstable modes, associated to elementary dynamical events such as bond 
breaking and reformation. On the contrary, the more fragile Lennard- Jones mixtures 
are characterized by a larger fraction of extended unstable modes, which lead to a 
more cooperative and heterogeneous dynamics. 
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1. Introduction 

Network-forming amorphous materials are of great interest for technological 
applications, as well as of fundamental importance for the theoretical understanding 
of the glass transition. At a microscopic scale, the structure of network glass-formers, in 
both the amorphous and highly viscous regime, is characterized by strong chemical 
ordering and atomic arrangements that usually form an open tetrahedral network. 
Upon cooling from high temperature, transport coefficients and structural relaxation 
times r of network liquids display a mild temperature dependence, often describable by 
the Arrhenius law r exp(E/T). Network glass-formers are thus "strong" in the 

Angell's classification scheme [1]. In contrast, other classes of glass-formers, including 
molecular, polymeric, and bulk metallic liquids, show super-Arrhenius temperature 
dependence of r, i.e., "fragile" behaviour. 

Ever since the introduction of the Angell's classification, the nature of the 
distinction between fragile and strong liquids has been highly debated. While the 
degree of fragility of a liquid correlates quantitatively with other macroscopic physical 
properties, the existence of qualitative differences between strong and fragile systems 
has been questioned. Evidence of a "fragile-to-strong" crossover in simulated network 
liquid [2] and numerical investigations of dynamic heterogeneities in model glass- 
formers [3] suggest that network liquids may just be an extreme case of the class of 
fragile systems [4]. In contrast, theoretical work on kinetically constrained models of 
glassy dynamics [5] indicates that strong behaviour may arise from the different nature 
of dynamical constraints. Moreover, the energy landscape description of supercooled 
liquids [6, 7] shows that the organization and connectivity of stationary points in the 
Potential Energy Surface (PES) may be qualitatively different in fragile and strong 
glass-formers. In this two latter scenarios, fragile and strong liquids would thus belong 
to different "universality classes" of glass-formers. 

Silica is often considered as a prototypical network glass-former. In recent years, 
several authors have studied structural and dynamical properties of this system through 
numerical simulations, employing both Molecular Dynamics (MD) and Monte Carlo 
techniques. One of the most realistic and popular models of silica available for molecular 
simulations is the BKS model introduced by Van Beest et al [8]. In this model, the 
interaction between Si and O atoms is described by a long-ranged Coulombic interaction, 
plus a short range repulsion of the Born-Mayer type. Various physical aspects of the 
supercooled and glassy regime of the BKS model have been analysed, including the phase 
diagram [9, 10], structural [11], dynamical [12, 13, 3, 4], and vibrational [14, 15, 16] 
properties. Investigations of the energy landscape of the BKS model have also been 
performed [17, 18, 19, 20]. Because of the long-ranged nature of the interactions, 
however, simulations using the BKS model are computationally demanding. Hence, 
development of simpler force-fields, capturing the basic features of network liquids, is 
highly desirable. Recently, in fact, simplified models for silica have been proposed, 
including short-ranged variants of the original BKS potential [21, 22] and primitive 
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models based on patchy interactions [23, 24]. Other models of tetrahedral network 
liquids (not directly related to silica) based on spherical, patchy interactions have also 
been studied recently [25] and in the past [26]. 

In this work, we present a new model of network glass-former, based on spherical, 
short-ranged potentials. Our model allows efficient simulations and can be tuned to 
reproduce some relevant properties of amorphous silica. It does not aim at a realistic 
description of liquid and amorphous silica, yet it captures to a good extent the essential 
physics of network glass-formers. Moreover, being able to describe both network and 
"simple" glass-forming liquids [27, 28] with similar efficiency via the same family of 
interactions, we can get an unusually detailed and systematic comparison between the 
microscopic origins of their structural relaxation. In particular, we trace back the 
distinct dynamic features of network glass-formers (e.g. strong behaviour, weak dynamic 
heterogeneity, bond breaking and reformation processes) to the properties of the PES, 
contrasting our findings with the case of the more fragile, close-packed Lennard- Jones 
(LJ) mixtures [27, 28]. Our results emphasize the role of the unstable modes of the 
PES, as a key to rationalize the different dynamic behaviours of glass- forming liquids. 

The paper is organized as follows: in section 2 we introduce our model of network 
glass-former; in section 3 we describe its structural and dynamical properties, while in 
section 4 we analyse the properties of the stationary points of the PES, focusing on the 
unstable modes. Finally, in section 5 we draw our conclusions. 



2. Model 



Our model of network glass-former, called NTW herein, is a binary mixture of classical 
particles interacting through the following forcefield 



cr r 
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u aa {r) = 4e aa — (1) 



u a p(r) = 4e a/3 



r 



a ^(3 (2) 



r 

where a, (3—1,2 are indexes of species. In the following, we will use an, en and 
y / m 1 cr^ 1 / 'en as reduced units of distance, energy and time respectively. Keeping an eye 
on silica, we identify large particles (species 1) with Si atoms and small particles (species 
2) with O atoms, and we fix the number concentrations at X\ = 0.33, x<i = 0.67. We also 
use the same mass ratio of O and Si atoms: m 2 /m 1 = 0.57. A smooth cut-off scheme 
is used to ensure continuity of u a/ 3(r) at r = 2.2<r Q( g up to the second derivative [29]. 
The size of the samples considered in this work is N = Ni + N2 = 500. The presence of 
finite size effects have been checked through simulations of larger systems (N = 2048, 
8000) and will be briefly discussed in section 3. We performed MD simulations in the 
NVE ensemble using quenching protocols and equilibration criteria similar to the ones 
of previous simulations of LJ mixtures [27, 28]. Equilibration and production runs 
were performed using Berendsen rescaling and velocity-Verlet algorithm, respectively. 
The time step St was varied between 0.001 (at high T) and 0.004 (at low T). The 
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Figure 1. Partial pair correlation functions g a p(r) for the NTW model (solid lines) 
and the BKS model for silica (open points). The thermodynamic state shown for BKS 
silica is p = 2.37g/A 3 ,T = 2750K and the one of NTW is p = 1.655, T = 0.39 in 
reduced units. The BKS data refer to the MD simulations by Horbach and Kob [30]. 



absence of major systematic aging effects was checked by comparing thermodynamic, 
structural, and dynamical properties in different parts of the production runs. At the 
lowest temperatures, simulations involved up to 3.5 x 10 7 and 7 x fO 7 steps for the 
equilibration and production runs, respectively. Thanks to the short range of the 
potentials and to the open local structure of the system, these long runs took a few 
days on a 3.4 GHz Xeon processor. 

To reproduce the open, tetrahedral local structure of network glasses, two main 
physical ingredients must enter in the forcefield of our model: highly non-additive 
interaction radii and strong attraction between unlike species. Building on previous 
experience [26], we determined the following optimal set of interaction parameters: 

012/011 = 0.49 022/011 = 0.85 
€l2 /e n = 6.00 e 22 /en = f.00 

To optimize the parameters above, we performed a series of preliminary simulations at 
reduced density p = p expt ~ 1.53, which corresponds to the density of amorphous silica 
in normal experimental conditions. The parameters were adjusted by requiring that the 
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0.020 




to [THz] 

Figure 2. Vibrational density of states (VDOS) obtained from local minima of 
the potential energy at p = 1.53, T = 0.30 for the NTW model (solid line), 
compared to experimental VDOS of amorphous silica (filled points, from [31]). The 
experimental data are convoluted with a correction function that alters their features 
only quantitatively [16]. 

ratio between the positions of the first peaks in gu{r) and gn(r) was equal to that of 
Si-0 and Si-Si interatomic distances of amorphous silica. We also checked that at low 
temperature the average coordination numbers, obtained from the integral of the radial 
distribution functions, were close to the ideal tetrahedral ones, i.e., Zx2 = 4, Z21 = 2. 

The physical units of our model can be be fixed to reproduce some relevant 
properties of experimental and realistic models of silica, such as the BKS model. We 
fixed the length scale axx so that the position of the first peak of gn(r) matched 
the mean Si-Si distance (i.e., 3.12A) of amorphous silica. In this way we obtained 
a xi = 2.84A. To fix the energy scale en, we compared the shape of the radial distribution 
functions g a 8{r) to the ones obtained by Horbach and Kob [30] for BKS silica at the 
state point p = 2.37g/A , T = 2750K (see figure 1). The corresponding density in 
reduced units is p = 1.655. A good overall agreement of the liquid structure is found 
around T = 0.39, from which we estimate en ~ 7000K. Finally, we fixed the time 
scale of our model by adjusting the mass scale mi so as to reproduce typical vibrational 
frequencies of amorphous silica. Following previous studies (see for instance [14, 15, 16]), 
we determined the vibrational density of states (VDOS) through diagonalization of 
the dynamical matrix calculated at local minima of the potential energy at p = 1.53, 
T = 0.30. The choice mi = 8.7 x 10 _23 g ~ 1.9msj yields reasonable agreement between 
the VDOS of our model and the experimental VDOS of amorphous silica [31] (see 
figure 2, which is further discussed in section 3). From the value of mi given above we 
obtain the time unit t Q = \fm^a\Jt\x = 2.0 x 10~ 13 s. 

3. Structure and dynamics 

In this section we further validate our model by analysing its structural and dynamical 
properties. Our simulations spanned a wide range of densities: 1.250 < p < 2.300. At 
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Figure 3. Partial structure factors Sn(k) (solid line), Si2(fc) (dashed line), and S22(k) 
(dotted line) at p = 1.655 and T = 0.29. 

higher density (p = 2.800) we found clear signs of crystallization of our samples, but 
we did not attempt to determine the crystallographic structure. At lower densities 
(p < 1.250) large voids are formed in the network structure, and liquid-gas phase 
separation might occur. In the following, we will mostly focus on the isochore p = 1.655, 
which corresponds to the density employed in several simulations of BKS silica, at 
temperatures in the range 0.29 < T < 1.50. 

3.1. Structure and vibrations 

The fact that the radial distribution functions of our model agree rather well with those 
of the more realistic BKS model (see figure 1) and the overall qualitative shape of the 
VDOS (see figure 2), already suggest that the NTW model should capture some relevant 
physical aspects of network glass-formers, at least for densities and temperatures where 
tetrahedral local ordering is more pronounced. In this section we study in more details 
the structural and vibrational properties of our model. 

In figure 3 we show the partial structure factors S a! 3(k) obtained at the lowest 
temperature attained in equilibrium conditions for p = 1.655. The pre-peak (also called 
first sharp diffraction peak) at k ~ 5.0 in Su(k) and S22{k) signals the formation of 
intermediate range order. This is a typical feature apparent at low temperature in 
network liquids. The positions of the pre-peak and main peak (k « 8.0) in Su(k) are in 
good agreement with those of SsiSi(k) in amorphous silica and simulated BKS silica [12]. 

Further insight into the structural properties of the NTW model is provided by the 
distribution f a i3 1 (9) of angles formed by a central particle of species (3 with neighbours 
of species a and 7, where a, (3, 7 = 1,2. Particles of species a and 7 are considered 
neighbours if their distance is less than the minimum of the radial distribution function 
at the corresponding T. The normalized angular distribution functions j '121 (0) and 
f2i2(&), shown in figure 4 for a few selected temperatures, reveal the typical features 
associated to local tetrahedral ordering. The broad peak in f 212(G), located around 
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Figure 4. Angular distribution functions / Q /3 7 (6') for T = 0.80 (dotted lines), 0.48 
(dash-dotted lines), 0.36 (dashed lines), 0.29 (solid lines) at density p = 1.655. 



6 = 108°, reflects the presence of slightly distorted tetrahedra centered around particles 
of species 1. The /i2i(#) shows a peak around 180°, which corresponds the links 
formed by particles of species 2 connecting adjacent tetrahedra. Note that the peak 
positions and the overall shape of these distribution functions change only mildly below 
T m 0.50. Thus, below this temperature, which we will identify in the next section as 
the onset temperature of slow-dynamics [32], the NTW model displays a strong degree 
of tetrahedral local ordering. 

The angular distribution functions fm(9) and j 222 (#) (see two lower plots in 
figure 4) provide information about the intermediate range order of our model. The 
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Figure 5. Fraction of particles with ideal tetrahedral coordination as a function of 
1/T: P(Zi2 = 4) (empty circles) and P(Z 2 i = 2) (filled circles) are shown for particles 
of species 1 and 2, respectively. The vertical dotted line marks the onset of the slow- 
dynamics regime. 



fm(9) displays a broad peak located at 9 « 105°, associated to distorted corner-sharing 
tetrahedra. The smaller peak around 60°, due to three-fold rings [33], decreases in 
height upon lowering the temperature. At higher density (p = 2.300, not shown here) 
this small peak increases in intensity (at fixed T), while the peak at 9 ~ 105° splits in two 
sub-peaks. Similar variations upon compression were found in the angular distribution 
functions of a more realistic model of silica [33]. Hence, we conclude that our simple 
model is able to capture some non-trivial structural features of network glass-formers. 

To highlight the formation of a nearly ideal tetrahedral network at low T, we show in 
figure 5 the T-dependence of the fraction of particles with ideal coordination numbers, 
i.e., P(Zi2 = 4) and P{Z 21 = 2) for particles of species 1 and 2, respectively. To 
identify neighbouring particles we used the same criterion as for the angular distribution 
functions discussed above. The fraction of ideally coordinated particles is already 
substantial around T « 0.5 [P(Z 12 = 4),P(Z 21 = 2) > 0.70] and approaches unity 
at low T. This provides further indication that for the density considered here the 
system is indeed in the optimal region of network formation [25]. 

A closer inspection of figure 2 shows that the VDOS of the NTW model reproduces 
all the qualitative features of the experimental VDOS of amorphous silica. The 
relative positions of the peaks in the VDOS of NTW match well enough those of the 
experimental VDOS. Note that the absence of a peak at small frequencies (u ~ 4 
THz) in the experimental data is due to insufficient experimental resolution [16, 34]. 
A careful comparison of simulated and experimental VDOS of silica can be found 
in [16]. Here we only recall that the VDOS of the BKS model is somewhat unrealistic 
at low and intermediate frequencies [16]. Similar deficiencies have also been found in 
recent modifications of the original BKS model employing short-ranged potentials [22] . 
Specifically, the distinct peaks at 12 and 24 THz, as well as the small peak around 18 
THz (D 2 line), are missing in the VDOS of BKS silica. Given the simplicity of the 




Figure 6. Intermediate scattering functions F s (k,t) (self part) at p = 1.655 for wave- 
vector k = 5.0 (top panel) and k = 8.0 (bottom panel). 



forcefield employed, the success of the NTW model in reproducing the main qualitative 
vibrational features of amorphous silica is rather remarkable. 

3.2. Relaxation dynamics 

Our first step in the description of the dynamical properties of the NTW model consists 
in the identification of the so-called "slow-dynamics regime" [32]. In this temperature 
regime, the dynamical properties of glass-forming liquids assume all their distinct 
features, including two-step relaxation, dynamic heterogeneities, etc. To detect the 
onset of slow-dynamics, we study the variation with temperature of the incoherent 
intermediate scattering functions 

, N a 

K(k,t) = — E< ex P^ • [»"*(*) " r >(°)]» ( 3 ) 

i=l 

where a = 1, 2 is an index of species. The t-dependence of F}(k, t) is shown in figure 6 
for temperatures in the range 0.29 < T < 1.50 at two different wave-numbers: k = 5.0, 
close to the pre-peak in the static structure factors (upper panel) and k = 8.0 (lower 
panel). Two-step relaxation develops around To ~ 0.50, which we take as the onset 
temperature of the slow-dynamics regime. Distinct damped oscillations are observed in 
F™(k,t) on the time scale of the early /3-relaxation, i.e., on approaching the plateau. In 
larger samples (not shown here), the amplitude of these oscillations is slightly smaller — a 
well-known finite-size effect in model network liquids [30, 35]. 
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Figure 7. Top panel: Angell plot for structural relaxation times r±(k = 5.0) (filled 
squares), r±(k = 8.0) (stars), and inverse diffusion coefficients l/Di (empty squares). 
All quantities refer to particles of species 1. The modified VFT fit for Ti(k = 5.0) 
[equation (4)] and the Arrhenius fit for l/Di arc also shown as solid lines. The 
dotted vertical line marks the onset of the slow dynamics regime. Lower panel: Ratios 
n(fc = 5.0)/r 2 (fc = 5.0) (filled squares) and D 2 /Di (empty squares) as a function of 
1/T. 



We now analyse the T- dependence of the structural relaxation times extracted from 
the intermediate scattering functions. Wave-number dependent relaxation times, r a (k), 
for species a are defined by the condition F"(k, T a (k)) — 1/e. In the Angell plot in Fig. 7 
we focus on the T-dependence of T\. We focus here on the case r = T\{k = 5.0). To fit 
the T-dependence of the relaxation times we use the following modified Vogel-Fulcher 
equation, previously employed in our study of LJ mixtures [27], 



r(T) 



Tooexp 
Tcoexp 



i 



K(T/T -l) 



T > T* 



where 



Tooexp 



Eoo/T* 



(4) 



(5) 



K(T*/T -1) 

Equation (4) describes the crossover from Arrhenius to Vogel-Fulcher T-dependence 
of r, ensuring continuity at T = T*. Its use in the case of a network glass- 
formers is justified by the observation that network liquids display a mild super- 
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Figure 8. Fragility index K obtained for NTW model and two representative LJ 
mixtures as a function of p (NTW, lower axis) and P (LJ, upper axis). 



Arrhenius behaviour around and slightly below Tq. As we can see from figure 7, 
equation (4) fits rather well t(T) over about 5 decades. Note that the degree of super- 
Arrhenius behaviour in r(T) is indeed rather modest and more visible at wave-numbers 
corresponding to the first sharp diffraction peak (k = 5.0). 

Also included in figure 7 are the partial diffusion coefficients D\(T) obtained from 
the usual Einstein relation. To describe the T-dependence of the diffusion coefficients, 
we simply used the Arrhenius law, D a = exp(E a /T). By fitting the data at 
low temperature (T < 0.4), we obtain activation energies E\ = 6.3 ~ 3.8eV and 
E 2 = 5.9 ~ 3.6eV, which are in reasonable agreement with those obtained in the 
case of BKS silica for silicon and oxygen atoms, respectively [30]. The difference in the 
diffusion coefficients between the two species is analysed in the lower panel of figure 7, 
where the ratio D1/D2 is shown as a function of 1/T. This ratio becomes ~ 2 at the 
lowest temperatures. A smaller separation of time scales is found when inspecting the 
ratio T2/V1. 

The fragility index K of the NTW model obtained from fits to equation (4) is shown 
in figure 8 as a function of p. The system is slightly stronger (smaller K) at densities close 
to the experimental density of silica. In this range of density (1.5 < p < 1.8), tetrahedral 
local order becomes nearly ideal at low T. Hence, our results provide support for the 
link between structure and dynamic behaviour in network liquids demonstrated in [24] 
for patchy colloidal particles. Interestingly, the fragility of the NTW model seems to 
increase outside the density range mentioned above both at high and low density. Further 
investigations at low density would be required to clarify the nature of this behaviour 
and its possible connection with a reversibility window [36, 37], whose existence for silica 
has been suggested by recent work [38]. 

The use of a common functional form to describe r = r(T) of NTW and L J models 
allows a direct comparison of their Angell's fragility. To this end, we also included 
in figure 8 the values of K obtained in [27] for the mixture of Kob and Andersen 
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(BMLJ [39]) and the mixture of Wahnstrom (WAHN [40]). Clearly, both LJ mixtures 
have larger fragility indexes. Furthermore, the fragility index of NTW, K = 0.09, 
obtained at p = 1.655 is lower by around a factor of 3 than the lowest values found for LJ 
mixtures (K = 0.24 for AMLJ-0.60 [27]). The fact that the fragility index was obtained 
at constant density for the NTW and constant pressure for the LJ mixtures does not 
affect substantially our conclusions. Moreover, even when considering the variation of 
K with p, the largest fragility index of the NTW model (K = 0.20 at p = 2.300) is 
comparable to the lowest ones found in LJ mixtures [27]. Thus, despite the presence 
of super-Arrhenius behaviour around the onset of slow dynamics, our network glass- 
former is stronger than all LJ mixtures studied in [27], a fact which fits naturally in the 
Angell classification scheme. Moreover, our analysis does not exclude the occurrence, at 
low temperatures, of a fragile-to-strong transition — a scenario which has recently found 
support on the basis of an energy landscape approach [20]. 



3.3. Dynamic heterogeneity 

Figure 8 shows that NTW, BMLJ, and WAHN may be considered as models of strong, 
intermediate, and fragile glass-formers, respectively. This offers the opportunity to 
investigate the main trends of variations of the dynamics in liquids with different 
fragility. In this section, we focus on the degree of dynamic heterogeneity of the above 
mentioned models. 

As a simple measure of the degree of heterogeneity of the dynamics we will use the 
non- Gaussian parameter 

3(r 4 (t)) 

Mi) = Tj 2mX 2 " 1 ( 6 ) 

which measures the deviation of the distribution of particles' displacements r(t) from 
a Gaussian distribution. Upon cooling the liquid below To, in fact, the distribution of 
particles' displacements deviates progressively from a Gaussian and the amplitude of 
a?2 increases. Within the late /5-relaxation time scale, the non-Gaussian parameter of 
typical glass-forming liquids displays a broad peak, whose the position t* and the height 
a* increase by decreasing temperature. The trends of variation of the maximum of the 
non-Gaussian parameter a\ have been found to follow qualitatively the behaviour of 
more refined dynamic indicators [3] , such as those obtained from four-point correlations 
functions [41, 42]. 

We computed the non-Gaussian parameter a^(t) in equation (6) separately for 
species 1 and 2. The results obtained for NTW, BMLJ, WAHN models along isochoric 
quenches are shown in Fig. 9 and 10 for species 1 and 2, respectively. The degree of 
dynamic heterogeneity, as measured from the height a\ of the peak, is least pronounced 
in the case of the NTW model close to the ideal density for tetrahedral network structure. 
This is consistent with the analysis of Vogel et al [3] , who found in fact that BKS silica 
had a lower degree of dynamic heterogeneity than other simple glass-formers, including 
the BMLJ model. Our results thus indicate a broad correlation between fragility and 
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the degree of dynamic heterogeneity in glass-forming liquids. In particular, within the 
slow-dynamics regime, both the local structure and the local dynamics appear more 
homogeneous in network than in close-packed glass-formers. 

3.4- Local rearrangements 

We now turn to a closer inspection of the nature of local rearrangements in our model 
network glass-former. In particular, we want to identify the structural modifications 
that accompany relaxation events. This is motivated by the current interest in 
investigating the link between structure and dynamics in glass- forming liquids [43, 44]. 
We will contrast the results for the NTW model to those previously obtained for LJ 
systems [45, 27, 28]. 

In a first attempt to characterize the local dynamics of the NTW model and to 
establish a connection with its local structural properties, we computed the "propensity 
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Figure 10. Same as Fig. 9 but for particles of species 2. 



of motion" of particles, according to the definition of Widmer-Cooper et al [43]. 
In this approach, time-dependent atomic displacements Ar(t), relative to a reference 
configuration, are averaged over several trajectories generated by independent initial 
sets of velocities ( "iso-configurational ensemble"). The resulting spatial distribution 
of average displacements (Ar(t))j C , where (. . .)j C denotes an average in the iso- 
configurational ensemble, is thus strictly associated to the initial configuration, and can 
be used, in principle, to identify the local structural features responsible for relaxation 
events. 

The spatial distribution of the particles with large propensity of motion for a 
representative configuration sampled at T = 0.31 is shown in figure 11 for t = 200 < t* 
(left panel) and for t = 1000 ~ t* (right panel). Mobile particles, depicted in figure 11 
as large dark spheres, are identified as the ones having the 30% largest propensities of 
motion among those of the same chemical species. The overall picture does not change 
upon small variation of the fraction of particles displayed. There is no substantial 
clustering of mobile particles on either time scales. Some clustering is observed at 
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(a) t = 200 



(b) t = 1000 «s t* 



Figure 11. Snapshots of particles having large propensity of motion (Ar(i)), c for 
t = 200 (left panel) and t = 1000 w t* (right panel). For both times, the 30% most 
mobile particles of either species are shown as large dark spheres, irrespectively of 
chemical species. The remaining particles are shown as small light spheres. 



t = 1000 but the size of the clusters remains rather modest (less than ~ 10 neighbouring 
particles). This is strikingly different from the results obtained in LJ systems [45] within 
the slow-dynamics regime. In L J systems, a pronounced clustering of particles with large 
propensity of motion has been observed for times on the order of the late /3-relaxation 
(t ~ t*). Our results confirm that the degree of dynamic heterogeneity of network 
liquids is much less pronounced than in close-packed LJ systems, and show that the 
origin of the weak dynamic heterogeneity observed within the a-relaxation time scale is 
essentially kinetic, rather than structural. 

The results above do not imply, however, that there is no link at all between local 
structure and dynamics in network liquids. Such a link is more subtle and requires 
a different type of investigation. The breaking and reformation of bonds involved 
in typical relaxation events [30] occurs, in fact, on a very short time scale and the 
averaging introduced by the iso-configurational ensemble washes out this information. 
To overcome this problem we calculated, following Ladadwa and Teichler [46], the 
instantaneous mobility of particles from a smoothed atomic trajectory 



where <f)(t',t) is a smoothing function normalized to 1. Rather than a Gaussian [46], 
we used a simple window smoothing function of length 2 At, which equals 1/(2 At) for 
t — At < t' < t — At and otherwise. From the smoothed trajectories, we computed 
the instantaneous atomic mobility [46] 




(7) 
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Figure 12. Instantaneous mobility m(t) in arbitrary units (dashed line) for two 
representative particles of species 1 (two upper panels) and species 2 (two lower panels) 
at T = 0.29. Also shown are the instantaneous coordination numbers Z\2(t) and ^i(t) 
for particles of species 1 and 2, respectively (solid lines). 



using At = 10. The time dependence of m(t) is shown in figure 12 for representative 
particles of species 1 and 2 at T = 0.29. At this low temperature, atomic mobilities 
show intermittent behaviour, with long periods of inactivity (vibrations) followed by 
displacements occurring on a very short time scale. To establish the connection with the 
changes in the local structure, we also plot, for the same time interval, the instantaneous 
coordination number Z(t). A clear correlation between intermittent dynamical events 
and bond breaking and reformation processes is observed. In particular, defective 
local environments (Z12 7^ 4 and Z21 7^ 2), either created instantaneously by thermal 
fluctuations or associated to long-lived defective configurations, are closely associated 
to dynamical events. Interestingly, this shows that in network liquids the link between 
structure and dynamics can be understood at a single-particle level. Such a link has 
been demonstrated in LJ systems only at a coarse-grained spatial level [47]. 




Figure 13. Animations of representative elementary dynamical events at T = 0.29 
(see http : //www-df t . ts . inf n . it/media/gp/1/ for the corresponding MPG files, 
realized with VMD [48]). White large spheres and small red spheres correspond to 
particles of species 1 and 2, respectively. The particles involved in the elementary 
dynamical events are surrounded by a yellow halo. Atomic positions have been 
averaged over a time window of 4.2 reduced time units (700 time steps) to remove 
thermal motion and help visualization. Left panel: four small particles (indicated 
by arrows in the figure) in the central-upper part of the figure perform a correlated, 
rotational motion in neighbouring tetrahedra (file moviel.mpg; size: 1.4 Mb). Right 
panel: the central particle (indicated by an arrow) explores its low density environment 
with a sequence of two jumps, associated to bond breaking and reformation (file: 
movie2.mpg; size: 1.7 Mb). 



Finally, we describe qualitatively the typical local rearrangements observed 
at low temperature in the NTW model. By inspection of animated atomic 
trajectories, we identified two typical relaxation processes, closely related to the ones 
occurring in BKS silica. See figure 13 and supplementary materials (available at 
http : / /www-df t . ts . inf n . it/media/gp/1/) for two representative events. A first class 
of rearrangements involves correlated rotations of tetrahedra formed by small particles 
around nearly immobile large particles. The overall process resembles the "rotational 
period" recently described by Heuer and coworkers [49], in which oxygens perform 
permutation of the tetrahedral positions around a fixed silicon atom. A second class 
of rearrangements is closely related to the ones described by Horbach and Kob [30]. 
One large particle jumps out from one of the faces of the tetrahedron surrounding it 
and attaches itself to an under-coordinated (Z 2 \ = 1) small particle. At the same 
time, a slight recoil movement of the small particles forming the involved tetrahedron is 
observed, together with the formation of a new dangling bond. Contrary to rotational 
rearrangements, which often involve a few neighbouring tetrahedra, this second class 
of elementary dynamical events is strongly localized around the involved tetrahedron. 
On longer time scales, however, sequences of independent events are also observed (see 
right panel of figure 13). 
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4. Stationary points and unstable modes 

Summarising our previous analysis, two key features characterize the dynamical 
behaviour of our model network liquid: strong behaviour in the Angell's classification 
scheme and a significant homogeneity of atomic displacements within the late (3- 
relaxation time scale. In this section, we rationalize these features in terms of the 
properties of the Potential Energy Surface (PES). In particular, we first provide an 
estimate of the average energy barriers in the PES of the NTW model, and then analyse 
the localization properties and real-space structure of the unstable modes associated to 
stationary points of the PES. 



4-1. Energy barriers 

The nature and distribution of barriers connecting stationary points of the PES have 
been long recognized as key aspects for understanding the dynamics of fragile and strong 
liquids [6, 7]. Recently, sophisticated analysis of transitions between meta-basins for 
model glass-forming liquids [50] have provided even further quantitative evidence of the 
importance of the PES. In this section, we employ a simple definition of average energy 
barriers [51, 28] to quantify the roughness of the potential energy surface of the NTW 
model. 

Our analysis of the PES is based on the procedure described in [28]. For each 
state point, we perform minimizations of the mean square total force W to locate the 
closest stationary points along the dynamical trajectory. Typically, between 100 and 400 
configurations per state point are considered as starting points for iy-minimizations. 
It is well known that iy-minimizations often locate points with a low value of W 
(W ~ 10~ 2 -T- 10~ 4 ) that are not true stationary points. These points, usually called 
quasi- saddles, contain nonetheless relevant information about the dynamics [52, 53]. In 
the following, we will include these points in our analysis, without further distinction 
between true stationary points and quasi-saddles. Having located the stationary points, 
we diagonalize the Hessian matrix of the potential energy surface and thus obtain a set 
of 3N eigenfrequencies u a and eigenvectors {e"}. The unstable eigenvectors (u^ < 0) 
are of particular interest for our discussion, because they are more directly related to 
the dynamical behaviour of the system [45]. As in previous work [28], we will report 
the imaginary branch of the frequency spectrum along the real negative axis. 

To estimate the average potential energy barriers we follow the definition given by 
Cavagna [51] 

E s = -— (8) 

where e s = e s (f u ) is the average energy of stationary points having a fraction of unstable 
modes f u = n u /3N. As in [28], we evaluate E s from the slope obtained by linear 
regression of e s vs. f u of individual stationary points sampled at temperature T. The 
procedure is illustrated in the upper panel of figure 14, where f u is shown as a function of 
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Figure 14. Upper plot: fraction of unstable modes f u as a function of the energy e s 
of individual stationary points sampled at the indicated temperatures. Also included 
are linear fits, from which the average energy barriers E s are obtained. Lower panel: 
average energy barriers E s as a function of the average energy of stationary points 
e s = e s (T) for various T. 

e s for stationary points sampled at selected temperatures (T < To). The energy barriers 
obtained for individual state points are collected in the bottom panel as a function of 
the average energy of stationary points e s = e s (T). 

From the comparison of the results above with those obtained in [28] for LJ 
mixtures, we draw two main conclusions, which highlight the peculiarity of the network 
liquids: (i) energy barriers in the NTW model are large compared to typical thermal 
energies already for T To, and (ii) their increase is very weak (less than 20%) with 
decreasing temperature below To- At least at a qualitative level, (ii) confirms our 
conjecture [28] that the fragility of a glass-forming liquid is related to the increase 
of average energy barriers E s upon cooling below To- It also supports the overall 
picture that the energy landscape of network liquids has a uniformly rough structure, 
with barriers whose average amplitude is nearly independent of the energy level. 
Organization of stationary points into meta-basin structures, while present even in 
network liquids [54], should be of much more limited extent than in the more fragile, 
close-packed glass-formers. 

4-2. Localization properties 

We now investigate the localization properties and the real-space structure of the 
unstable eigenvectors of the stationary points sampled in the slow-dynamics regime. 
In particular, we aim at explaining the more "homogeneous" character of atomic 
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Figure 15. Average gyration radius L c of modes of frequency w calculated at p = 1.655 
and T = 0.31 using different procedures (see text for definitions): straightforward 
calculation using equation (9) (dotted line) , method (i) (dashed line) , and method (ii) 
(solid line). 



displacements observed in NTW, compared to the more fragile LJ mixtures. 

One possible measure of the degree of mode localization, used in a number of 
previous investigations [55, 56, 57], is provided by the gyration radius 

N 



where r c = X]i=i r *il e fl * s the "centre of mass" of the mode. Extended modes should 
have L c « 1.0, while L c < 1.0 for localized modes. It turns out, however, that the 
definition in equation (9) is inappropriate for systems with periodic boundary conditions. 
The centre of mass of the mode, in fact, is not well defined in a periodic system, and the 
value of L c thus depends on the choice of the origin of the frame of coordinates. As it will 
be clear in the following, this shortcoming is particularly evident in the case of strongly 
localized modes. To overcome this problem, we employ two alternative definitions of 
the gyration radius, obtained by redefining the origin of the system coordinates: (i) the 
position of the particle that has the largest displacement on mode a is used as the origin 
of the system coordinates for the calculation of L°; (ii) the gyration radius is determined 
by minimization of L° over all possible origins of the system coordinates chosen on a 
grid of points subdividing the cubic cell. 

In figure 15 we show the average gyration radius L c of modes with frequency uj 
obtained for NTW at T = 0.31, using the original definition and the two alternatives 
(i) and (ii). The original definition substantially overestimates the extension of strongly 
localized modes, while discrepancies are somehow less pronounced for extended modes. 
For extended modes, smaller discrepancies are apparent between methods (i) and (ii). 
In the following, we will employ definition (ii). Only minor quantitative differences in 
the following analysis appear when using definition (i). 

We now focus on the unstable modes, which have been found to contain direct 
information on the dynamics of a glass-forming liquid [45]. In figure 16 we show the 




(9) 



i=l 
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Figure 16. Distribution of L c for unstable modes sampled in the slow-dynamics 
regime: in NTW at p = 1.655, T = 0.31 (upper panel), in BMLJ at p = 1.2, T = 0.45 
(middle panel), in WAHN at p = 1.297, T = 0.54 (lower panel). 



distribution of L c for the NTW model at T = 0.29. For comparison, we also show 
analogous distributions obtained for BMLJ and WAHN at temperatures deep in the 
slow-dynamics regime. The distribution of L c for NTW is bimodal, with an excess 
of localized unstable modes having L c ~ 0.2. A similar, yet much less pronounced, 
excess peak at low L c is observed in the distribution for BMLJ, while no such feature is 
found for the very fragile WAHN model. Justified by the fact that the minimum of the 
distribution of L c for NTW is located around 0.5, we define a mode localized (extended) 
if L c is smaller (larger) than 0.5. The precise location of this cutoff is irrelevant for the 
discussion below. 

Information about the relevance of localized unstable modes in different glass- 
formers is conveniently represented by the plot in figure 17, where the fraction of 
localized unstable modes f u \ is shown as a function of f u . The NTW model has a 
substantial fraction of localized unstable modes, independent of the energy level in 
the PES. /„/ only weakly increases by lowering f u , i.e., as the systems explores lower 
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Figure 17. Fraction of localized unstable modes f u i in stationary points as a function 
of the total fraction unstable modes f u in NTW at p = 1.655 (squares), BMLJ at 
p = 1.2 (circles), and WAHN at p = 1.297 (triangles). 



regions of the energy landscape. In contrast, in the very fragile WAHN f u i is very 
small and increases only on approaching the bottom of the energy landscape. The 
BMLJ displays an intermediate trend, since a larger fraction of localized unstable modes 
is present. These localized modes of the BMLJ mixture correspond to eigenvectors 
where a small particle and few other neighbours have large displacements (see figure 13 
in [28]). It is remarkable that the trend observed in the localization of unstable modes 
follows the different dynamic character (strong, intermediate, fragile) of the models 
studied (see also [58]). Moreover, our results provide a simple explanation of the weak 
degree of dynamic heterogeneity in network liquids in terms of an excess of localized, 
uncooperative unstable modes, which are absent, or at least rarer, in the more fragile 
LJ mixtures. 

Finally, we describe the nature of the rearrangements associated to localized and 
extended unstable modes of the NTW model. The typical extensions of localized and 
extended unstable modes are depicted in figure 18. From inspection of the real-space 
structure of the atomic displacements we found that the unstable modes reproduce the 
two classes of local relaxation processes described in section 3. Extended unstable modes 
usually involve coupled rotations of tetrahedra and have a marked collective nature. 
These modes are thus good candidates for explaining the rotational motions described 
in section 3. Similar unstable modes have been found in the unstable branch of the 
instantaneous normal mode spectrum of BKS silica [19]. On the other hand, localized 
unstable modes correspond rather well to the second class of local rearrangements 
described in section 3. The mode depicted in left panel of figure 18 shows the passage of a 
large particle, initially at the centre of a tetrahedron, through the face of the tetrahedron, 
while a neighbouring under-coordinated small particle moves in the opposite direction 
to create a bond. In the intermediate stage along the reaction coordinate of the 
mode the large particle is 5-fold coordinated. Our results indicate that both types of 
elementary dynamical events should be taken into account for a complete description of 
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Figure 18. Typical extension of localized (a) and extended (b) unstable modes of a 
stationary point sampled at T = 0.29. Only particles with displacement \ef\ larger 
than 0.04 are shown, and the eigenvectors are scaled logarithmically. Particles of 
species 1 and 2 are shown as white large spheres and small dark spheres, respectively. 



the dynamics in network liquids. Approaches that focus only on soft stable modes [59, 60] 
may not be able to capture the localized nature of the dynamics in network liquids. In 
these systems, in fact, the low frequency portion of the VDOS encompasses collective 
modes that typically involve coupled rotations of tetrahedra [61, 14]. 

5. Conclusions 

In this work we have extended our previous analysis [27, 28] on the glass transition 
of fragile Lennard- Jones mixtures by introducing a new model of tetrahedral network 
glass-former based on short-ranged, spherical interactions. Remarkably, these simple 
models of liquids, all based on pair potentials of the Lennard- Jones type, are able to 
reproduce qualitatively a wide spectrum of dynamic behaviours, thus allowing extensive 
and detailed investigations of the glass transition phenomenon. 

Notwithstanding the problem of crystallization, which may occur in binary mixtures 
during longer simulations [62, 63] but is not observed in our samples, we have found that 
the fragile vs. strong behaviour of our models can be clearly identified and rationalized 
even at relatively high T, but below the onset temperature of slow dynamics. Using 
an appropriate parameterization of the T-dependence of relaxation times, we have 
found that the model network glass-former is stronger at all studied densities than 
all previously investigated LJ mixtures. Our results also confirm that the degree of 
dynamic heterogeneity is less pronounced in network than in close-packed glass-formers. 

An important aspect of the glass transition concerns the nature of atomic 
rearrangements occurring within the a-relaxation time [64, 65]. The relation between 
dynamical events and the nature of the local structure is of particular interest [43, 44]. 
The analysis of the propensity of motion of particles within the late /3-relaxation time 
scale, combined with a comparative study of the non-Gaussian parameter in different 
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systems, has revealed a substantial homogeneity of atomic mobility in the model network 
glass-former. However, contrary to the case of LJ mixtures, it is possible to establish a 
relation between local structure and dynamics at the single-particle level by considering 
individual atomic trajectories. Periods of high mobility are in fact clearly associated to 
sequences of bond breaking and reformation, i.e. variations in the local structure. 

The features described above and the variation of dynamic behaviour in systems 
with different fragility can be rationalized well in terms of the features of the Potential 
Energy Surface. We have focused on the properties of the unstable modes of saddles 
and quasi-saddles sampled within the slow-dynamics regime. The amplitude of the 
average energy barriers E s in the model network glass-former is always larger than 
typical thermal energies below Tq and depends very mildly on the energy level. This 
contrasts the findings in the more fragile LJ mixtures, where E s rapidly increases upon 
entering in the slow-dynamics regime [28]. The localization of the unstable modes offers 
direct insight into the elementary dynamic events leading to relaxation. In general, as 
the system explores lower and lower regions of the energy landscape, the unstable modes 
soften and retain a cooperative character. In the NTW model, there is also a significant 
fraction of localized unstable modes that persists in the whole slow-dynamics regime. 
These localized modes typically describe bond breaking and reformation, i.e. elementary 
rearrangements that characterize the dynamics of the model. On the contrary, close- 
packed fragile liquids have a large fraction of extended unstable modes, which soften 
and tend to localize only on approaching the bottom of the landscape. As a result, the 
dynamics in the latter systems is inherently more cooperative than in network liquids. 
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